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MR  Elastography  System  for  Breast  Cancer  Detection 


1  INTRODUCTION 

Breast  cancer  is  the  most  commonly  diagnosed  cancer  in  women.  Early  diagnosis,  which 
is  critical  for  favorable  clinical  outcomes,  is  complicated  by  the  difficulty  of  detecting  small 
tumors.  One  physical  property  that  clearly  distinguishes  healthy  from  cancerous  tissue  is 
mechanical  stiffness  (hardness).  For  this  reason,  palpation  has  long  been  used  for  early 
screening.  Researchers  have  attempted  to  combine  external  mechanical  stimulation  and 
Magnetic  Resonance  Imaging  (MRI)  to  quantitatively  measure  the  Young’s  modulus  of  tissue 
throughout  the  breast  (and  the  prostate  as  well).  This  technique,  Magnetic  Resonance 
Elastography  (MRE)  has  been  called  “palpation  at  a  distance.”  One  of  the  most  challenging 
technical  issues  of  MRE  is  the  solution  of  the  “inverse  problem,”  i.e.,  robustly  and  accurately 
estimating  the  distribution  of  Young’s  modulus  from  MRI-measured  tissue  displacement  data  in 
a  practicable  computing  time.  Addressing  this  problem  was  the  focus  of  this  project. 

Based  on  understanding  developed  over  the  course  of  this  project,  as  well  as  recent 
technical  advances  made  both  here  and  at  other  institutions,  the  emphasis  placed  on  different 
aspects  of  the  scope  of  work  changed  somewhat  from  what  was  originally  anticipated.  The 
following  summarizes  how  we  addressed  the  tasks  outlined  originally  in  the  proposal: 

Task  1:  Develop  a  Finite  Element  Model  for  the  Behavior  of  Breast  Tissue 

Considerably  greater  work  was  devoted  to  this  task  than  originally  anticipated.  Primarily, 
this  occurred  because  we  determined  that  much  of  the  work  published  in  the  literature,  which  we 
initially  used  as  a  basis  for  own  study,  assumed  incorrect  (too  small)  values  for  the  Poisson’s 
ratio.  A  reduced  value  improves  numerical  performance  but  has  the  effect  of  introducing 
artifacts  in  the  predicted  response.  For  this  reason,  two  sets  of  finite  element  codes  were 
developed  over  the  course  of  the  project:  a  conventional  formulation  based  on  compressible,  2-D 
continuum  mechanics,  and  an  innovative  model  of  incompressible  behavior.  The  latter  is 
considered  more  accurate  and  also  potentially  offers  substantial  gains  in  computational 
efficiency.  This  work  is  detailed  in  Sections  2.2,  2.3,  and  2.4. 

Task  2:  Develop  Inverse  Model 

This  work  proceeded  essentially  as  originally  planned,  and  the  effort  culminated  in  an 
efficient  means  for  solving  the  inverse  problem  using  an  innovative  “adjoint”  formulation.  This 
work  is  described  in  Section  2.5. 

Task  3:  Develop  Hardware  for  Mechanically  Exciting  Tissue 

Several  groups  were  actively  researching  MRE  during  the  time  period  of  the  project.  We 
determined,  based  on  reviews  of  the  literature,  that  several  approaches  were  being  used  to  excite 
tissue  in  a  manner  compatible  with  the  requirements  of  the  MRI  environment.  In  particular, 
piezoelectric  actuators  and  modulated  radiation  pressure  from  high-powered  ultrasound 
machines  appear  to  have  been  quite  successful.  Based  on  the  success  of  these  other  efforts,  as 
well  as  our  findings  that  more  significant  needs  exist  on  the  analytical  front,  we  concentrated  our 
efforts  on  Tasks  1  and  2. 
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2  BODY 

2.1  Background 

Breast  cancer  is  the  most  commonly  diagnosed  cancer  in  women,  with  a  current  mortality 
rate  of  approximately  23  per  100,000.  Early  diagnosis,  which  is  critical  for  favorable  clinical 
outcomes,  is  complicated  by  the  difficulties  associated  with  detecting  small  tumors.  While 
mammography  is  very  useful,  the  X-ray  attenuation  properties  of  tumors  and  normal  tissues  are 
similar,  often  making  it  difficult  to  distinguish  small  tumors  from  healthy  tissue.  The  same 
problems  afflict  ultrasound,  since  the  acoustic  impedance  of  tumors  is  only  slightly  different 
from  that  of  normal  tissue. 

One  physical  characteristic  that  does  clearly  distinguish  healthy  from  cancerous  tissue  is 
mechanical  stiffness  (hardness).  The  Young’s  modulus  of  breast  tumors  can  be  one  to  two 
orders  of  magnitude  larger  than  that  of  normal  tissue  (Muthupillai  and  Ehman,  1996).  For  this 
reason,  palpation  has  long  been  used  by  both  physicians  and  patients  for  early  screening. 
However,  manual  palpation  is  not  very  effective  for  tumors  lying  deep  within  the  breast  and  is 
not  quantitative.  For  these  reasons,  numerous  efforts  have  been  under  way  over  the  past  decade 
to  combine  external  mechanical  stimulation  and  various  imaging  techniques  to  quantitatively 
measure  the  Young’s  modulus  of  tissue  throughout  the  breast  (and  the  prostate,  as  well).  This 
has  been  termed  “palpation  at  a  distance.”  Because  of  its  inherent  sensitivity,  most  of  the  current 
research  is  focused  on  exploiting  the  capabilities  of  magnetic  resonance  imaging  (MRI),  and  the 
overall  process  is  termed  magnetic  resonance  elastography  (MRE). 

Probably  the  most  challenging  technical  aspect  of  MRE  is  the  solution  of  the  inverse 
problem,  i.e.,  quantitatively  inferring  Young’s  modulus  from  MRI-measured  tissue  displacement 
data  (Trahey,  2001).  Addressing  this  problem  was  the  focus  of  the  project. 

2.1.1  Previous  Work  on  MRE 

A  number  of  researchers  have  investigated  use  of  displacement  images  to  infer  the 
Young’s  modulus  of  tissue.  Early  work  focused  on  utilizing  ultrasound  images,  but  much 
current  research  employs  MRI,  which  can  accurately  measure  much  smaller  displacements  (less 
than  one  micron)  with  higher  resolution  than  ultrasound. 

MRE  techniques  can  be  broadly  distinguished  by  whether  they  employ  static  or  dynamic 
displacement  of  the  tissue  surface.  Static  techniques  measure  the  change  in  tissue  displacement 
resulting  from  a  quasi-static  push  on  the  outside  of  the  breast,  or  a  sequence  of  pushes.  This 
method  is,  in  principle,  simpler  than  dynamic  methods,  but  certain  mathematical  difficulties  arise 
from  the  use  of  static  displacements.  These  difficulties  appear  to  have  led  to  the  virtual 
abandonment  of  this  technique  (Raghavan  and  Yagle,  1994;  Chenevert  et  al.,  1998; 
Skovoroda  et  al.,  1995). 

Dynamic  MRE  usually  involves  the  generation  of  oscillatory  motions  at  the  surface  of  the 
breast  to  induce  shear  waves  in  the  tissue.  Because  shear  waves  are  strongly  attenuated, 
researchers  have  also  investigated  the  use  of  ultrasound  to  induce  shear  waves  in  deep  tumors  via 
radiation  pressure  (Wu  et  al.,  2000).  In  either  case,  the  motion  of  the  interior  of  the  breast  (or 
phantom)  is  then  measured  using  MRI,  and  the  elastic  modulus  is  inferred  from  these 
measurements. 
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MRE  requires  the  use  of  relatively  low  excitation  frequencies  since  motion  at  frequencies 
higher  than  about  1  kHz  cannot  be  resolved  by  MRI  due  to  limitations  in  the  maximum 
modulation  rate  of  the  gradient  field.  Shear  waves  are  used,  because  even  at  these  low 
frequencies  they  have  a  sufficiently  short  wavelength  to  allow  the  resolution  of  small  features 
such  as  tumors.  Longitudinal  waves  have  a  much  faster  sound  speed,  a  correspondingly  larger 
wavelength,  and  thus  little  resolving  power  at  low  frequency. 

Three  techniques  have  been  used  to  infer  Young’s  modulus  from  shear-wave-induced 
displacement  data.  Early  research  efforts  employed  direct  measurements  of  the  wavelength  of 
shear  waves  at  various  points  in  a  tissue  phantom  to  determine  shear  wave  speed.  From  shear 
wave  speed,  the  modulus  can  be  readily  determined  (Muthupillai  et  al.,  1995; 
Manduca  et  al.,  1996).  A  difficulty  of  this  direct  technique  is  that  in  real  tissue,  acoustic  waves 
scatter  off  tissue  inhomogeneities  and  the  scattered  waves  interfere,  generally  leading  to  a  very 
complex  pattern  of  displacements  (van  Houten  et  al.,  1999).  Such  a  complex  pattern  does  not 
lend  itself  to  direct  measurements  of  wavelength.  To  avoid  this  problem,  researchers  have  turned 
to  physics-based  modeling  of  tissue  motion.  In  one  implementation  of  this  method,  the  equations 
of  motion  of  the  tissue  are  defined  in  terms  of  the  known  oscillatory  boundary  conditions  and  a 
set  of  assumed  tissue  properties  (i.e..  Young’s  modulus).  These  equations  are  then  solved, 
usually  with  finite  element  techniques,  and  then  the  calculated  displacement  patterns  are 
compared  to  those  measured  with  MRI.  The  assumed  distribution  of  Young’s  modulus  in  the 
model  is  then  systematically  altered  and  the  process  repeated  until  the  predicted  and  measured 
displacement  patterns  match  to  within  some  tolerance.  In  another  implementation,  one  derives 
from  the  equations  of  motion,  a  partial  differential  equation  for  the  Young’s  modulus  in  terms  of 
input  tissue  displacements  (Bishop  et  al.,  2000). 

One  of  the  most  active  groups  pursuing  physics-based  approaches  is  at  Dartmouth 
College’s  Medical  School  and  Thayer  School  of  Engineering.  They  have  reported  substantial 
success,  but  two  key  problems  stand  out: 

1.  The  time  required  to  converge  the  inverse  calculations  can  be  impracticably  long, 
even  for  a  relatively  coarse  mesh  on  a  high-performance  workstation.  While  this  has 
led  to  the  development  of  so-called  subzone  techniques  and  utilization  of  parallel 
processing  to  reduce  execution  times,  this  continues  to  represent  a  significant 
practical  problem  for  using  MRE  in  a  clinical  setting  (Van  Houten  et  al.,  1999;  Van 
Houten,  et  al.,  2000a,  Van  Houten  et  al.,  2001). 

2.  While  good  results  have  been  reported  in  numerical  simulations,  the  results  in  real 
tissue  and  in  phantoms  have  been  decidedly  mixed.  This  has  lead  to  recent  efforts  to 
model  transient  effects  (rather  than  assuming  steady-state),  motion  in  three 
dimensions  (rather  than  the  plane  strain  approximation  used  previously),  and 
speculation  that  viscoelastic  behavior  may  need  to  be  modeled  (Van 
Houten  et  al.,  2000b).  Use  of  these  more  sophisticated  models  will  further  exacerbate 
the  run-time  problem. 

2.2  General  Aspects  of  Solving  the  Tissue  Equations  of  Motion 

Finite  element  techniques  based  on  physics-based  equations  of  motion  can  be  developed 
to  either  solve  for  the  tissue  displacements  given  an  assumed  set  of  Young’s  moduli,  or  to  solve 
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directly  for  the  spatially-dependent  stiffness  by  employing  a  complete  set  of  displacement  data. 
The  latter  approach  is  potentially  faster  but  tends  to  be  more  sensitive  to  measurement  errors,  and 
our  efforts  on  this  project  focused  on  the  former  scheme. 

To  develop  a  solution  for  the  inverse  problem  using  this  approach,  it  is  therefore  first 
necessary  to  develop  a  numerical  model  for  the  tissue  displacements  given  an  assumed  set  of 
tissue  properties,  the  so-called  “forward  problem.”  While  a  commercial  product  such  as  ANSYS 
can  (and  was)  used  for  preliminary  and  scoping  analyses  as  described  later,  ultimately  a 
customized  computer  code  is  needed  so  that  the  displacement  solution  can  be  integrated  with  the 
solution  of  the  inverse  problem. 

2.2.1  General  Form  of  the  Tissue  Equations  of  Motion 

For  dynamic  simulation  of  vibrations  in  tissue,  finite  element  codes  can  be  configured  to 
calculate  the  steady-state  amplitude  of  a  damped,  linear  elastic  solid  in  either  three  dimensions  or 
in  plane  strain.  Typically  the  boundary  conditions  are  specified  as  a  harmonically-varying  shear 
displacement  along  one  or  more  of  the  edges.  In  terms  of  developing  a  finite-element  model,  the 
equations  of  motion  of  the  nodes  of  the  model  can  be  written: 

KMu  +  0KM  —  +MM^-  =  F  (1) 

dt  dt2 


where 

u  =  time-dependent  displacement  vector  of  the  nodes 
KM  =  stiffness  matrix  (relates  displacements  to  stresses) 

MM  =  mass  matrix  (defines  inertia) 

P  =  Rayleigh  damping  coefficient  that  defines  the  damping  ratio  in  terms  of  the  stiffness  matrix 
F  =  time-dependent  applied  external  forces 


It  can  be  shown  that  if  the  desired  damping  ratio  is  R,  then  for  an  angular  excitation  frequency  co 
one  should  use  (Smith  and  Griffiths,  1998): 


m 


It  should  be  noted  that  it  is  also  possible  to  define  the  damping  in  terms  of  the  mass  matrix,  using 
a  slightly  different  form  for  the  coefficient. 

For  the  case  where  we  assume  that  a  steady-state  harmonic  excitation  at  angular  rate  to  is 
employed,  we  define  the  time-dependent  forces  and  motion  in  terms  of  their  amplitudes  via: 
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Substituting  Equations  (3)  and  (4)  into  Equation  (1),  we  obtain: 
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Equation  (5)  is  a  time-independent  matrix  equation  for  the  displacement  amplitudes  in  terms  of 
the  amplitudes  of  the  applied  forces  and  an  assumed  set  of  tissue  properties. 

The  details  of  Equation  (5)  will  differ,  depending  on  the  precise  set  of  assumptions  that 
are  made:  compressible  versus  incompressible,  two  versus  three  dimensions,  etc.  Before 
considering  these  details,  we  first  consider  general  aspects  of  solving  these  equations  that  are 
common  to  all  of  these  approaches. 

2.2.2  Numerical  Solution  of  the  Equations  of  Motion 

It  is  convenient  to  write  Equation  (5)  in  the  standard  matrix  equation  form: 

A **k=bt  (6> 

where  A  is-the  (complex)  matrix  involving  the  terms  in  parenthesis  in  Equation  (5),  b  represents 
the  (real)  amplitude  of  the  forces  supplied  at  the  surface,  and  x  is  the  solution  vector  which  is 
therefore  also  complex.  The  Einstein  summation  convention  is  used,  i.e.,  a  summation  is  implied 
on  the  repeated  index  k.  Physically,  the  fact  that  x  is  complex  arises  from  differences  in  phase 
between  the  displacements  at  different  points  within  the  solid.  As  evident  from  Equation  (5),  the 
variations  in  phase  are  caused  by  the  effects  of  damping.  Careful  inspection  of  the  terms  in 
Equation  (5)  will  also  reveal  that  the  matrix  A  in  Equation  (6)  is  symmetric. 

Actually,  one  ordinarily  does  not  know  the  forces  applied  at  the  surface  of  the  breast  but 
rather  the  magnitude  of  the  forced  displacements.  For  this  reason,  Equation  (6)  is  rewritten  in 
practice  so  that  the  surface  nodes  are  forced  to  have  a  known  vibration  amplitude.  One  way  to 
do  this  is  to  replace  the  associated  matrix  elements  in  Equation  (6)  as  follows: 

Aik  -  8ik  (for  all  nodal  displacements  i  subjected  to  external  forces)  (7) 

bj  =  amplitude  of  excitation  of  ith  degree-of-freedom  (8) 

The  terms  in  b  corresponding  to  non-extemally  forced  nodes  are  zero. 

However,  it  can  be  advantageous  to  use  an  even  simpler  scheme  in  which  the  on-diagonal 
element  An  and  the  corresponding  displacement  bj  of  the  equation  corresponding  to  a  known 
displacement  m,  are  both  multiplied  by  a  very  large  number,  so  that  the  external  forcing 
overwhelms  the  effects  of  internal  forces.  This  technique  preserves  the  inherent  symmetry  of  the 
matrix  A. 

A  crucially  important  aspect  of  the  inverse  problem  is  obtaining  an  efficient  solution  of 
the  matrix  Equation  (5)  or  (6).  It  should  be  recognized  from  the  outset  that  efficiency  is  more 
than  simply  a  convenience,  since  an  iterative  solution  of  the  inverse  problem  will  usually  involve 
a  very  large  number  of  solutions  of  the  forward  problem.  While  numerous  standard  packages 
exist  for  efficiently  solving  sparse,  complex,  Hermitian  matrices,  achieving  efficiency  is 
complicated  in  our  case  because  the  large  matrix  A  is  (at  best)  symmetric,  and  therefore  is  not 
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Hermitian  (i.e.,  A*  A+).  However,  a  direct  solution  scheme  is  available  for  equations  of  the 
form  of  (6)  that  takes  advantage  of  both  the  sparseness  as  well  as  the  symmetry  of  A 
(Harwell,  undated).  This  is  a  specialization  of  the  typical  LU-decomposition  scheme  ordinarily 
used  to  solve  real  matrix  equations  of  the  form  of  Equation  (6).  Aside  from  raw  efficiency,  a  key 
advantage  of  such  schemes  is  that  once  the  matrix  A  has  been  decomposed,  e.g.,  factored  into 
lower-  and  upper-triangular  matrices  L  and  U,  solutions  for  different  forcing  vectors  b  can  be 
obtained  with  virtually  no  additional  effort.  As  demonstrated  below,  this  is  very  useful  for 
calculating  the  gradient  of  functions  of  the  displacements  with  respect  to  the  input  parameters. 
In  any  event,  the  forward  problem  solution  of  the  symmetric  matrix  is  obtained  in  about  one- 
eighth  the  time  required  for  the  solution  of  an  asymmetric  matrix  of  the  same  size. 

As  part  of  this  project,  so-called  “mesh-free”  techniques  were  also  investigated.  When 
coupled  to  iterative  matrix  solution  techniques  (Smith  and  Griffiths,  1998),  the  advantage  of  this 
approach  is  that  the  complete  stiffness  and  mass  matrices  MM  and  KM  need  never  be  assembled 
into  memory.  This  allows  very  large  problems  to  be  solved  even  when  memory  is  relatively 
limited.  Good  results  were  obtained  by  solving  Equation  (6)  using  a  variant  of  the 
preconditioned  conjugate  gradient  technique  (PCG),  namely  the  biconjugate  gradient  squared 
iterative  method  (Joubert  et  al.,  undated).  For  example,  the  solution  of  a  20,000  element  grid 
with  20,000  nodes  was  obtained  in  less  than  about  30  seconds  on  a  1  GHz;  personal  computer. 
However,  when  sufficient  memory  is  available,  the  iterative  techniques  cannot  compete  with  the 
direct  approach,  so  in  practice  it  was  found  to  be  better  to  use  the  latter. 

2.3  The  Equations  of  Motion  of  a  Compressible  Medium 

In  this  section,  we  apply  the  general  framework  developed  in  Section  2.2  to  the  specific 
case  of  a  linear,  compressible  medium  in  the  plane  strain  approximation.  We  will  then 
investigate  some  problems  associated  with  the  compressibility  of  the  medium. 

2.3.1  Development  of  a  Finite  Element  Code 

To  support  numerical  experiments,  a  solution  of  Equation  (5)  was  developed  for  the 
compressible  equations  of  motion  in  the  plane  strain  approximation  using  standard  finite  element 
techniques  (Smith  and  Griffiths,  1998).  The  stiffness  and  mass  matrices  were  symmetric, 
allowing  the  efficient  Harwell  matrix  solution  package  to  be  used  for  the  solution. 

An  exact  solution  of  the  displacements  is  available  for  the  special  case  of  a  homogenous 
rectangular  solid,  fixed  on  all  but  one  edge,  and  subjected  to  a  transverse  excitation  whose 
amplitude  varies  sinusoidally  in  space  across  the  unfixed  edge  (Gao,  1995).  This  solution  was 
compared  to  the  predictions  of  the  Creare  finite  element  code.  The  latter  involved  a  10  cm 
x  10  cm  square,  modeled  using  1  mm  x  1  mm  3-node  triangular  elements.  Figure  1  shows  the 
displacement  amplitude  parallel  to  the  direction  of  excitation  down  the  midline  of  the  solid.  The 
code  results  agree  reasonably  well  with  the  exact  solution.  It  is  noteworthy  that  the  value  of  the 
damping  ratio  used  in  this  example  problem  (R  =  0.1)  is  not  quite  sufficient  to  prevent  some 
standing  waves  from  developing  along  the  side  of  the  rectangle  that  is  opposite  the  driven  edge. 

Of  course,  the  accuracy  of  the  solution  is  controlled  by  the  fineness  of  the  grid  and  the 
fidelity  afforded  by  the  various  types  of  elements.  Figure  1  utilized  a  fairly  coarse  mesh,  and 
more  accurate  solutions  can  be  obtained  by  using  more  elements,  quadrilateral  elements  instead 
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of  3-node  triangles,  or  triangles  with  more  than  3  nodes.  Figure  2  shows  the  behavior  of  a 
snapshot  of  the  displacements  at  a  particular  time,  i.e.,  the  real  part  of  the  complex  displacements 
(not  their  total  amplitude  as  shown  in  Figure  1).  The  cases  with  a  finer  mesh  (i.e.,  more  nodes) 
show  much  better  agreement  with  the  exact  solution.  This  experience  illustrates  an  important 
point:  the  discretization  of  the  problem  domain  must  be  adequate  to  resolve  tissue  motion  over 
the  wavelengths  that  will  develop  for  the  given  excitation  frequency.  Since  the  solution  time  will 
scale  as  the  cube  of  the  number  of  grid  points,  achieving  adequate  accuracy  in  a  practicable  time 
places  a  real  burden  on  the  analyst.  If  the  analysis  is  extended  to  three  dimensions,  the 
computational  burden  can  easily  become  very  severe  indeed. 

2.3.2  The  Effect  of  Poisson’s  Ratio  on  Sound  Speed 

The  example  calculations  with  the  Creare  finite  element  code  mentioned  above  used 
prototypical  values  of  the  Young’s  modulus  for  normal  breast  tissue  and  for  tumors.  However, 
the  value  of  Poisson’s  ratio  v  in  these  initial  calculations  was  set  to  0.3,  a  typical  value  for 
structural  (non-biological)  materials.  In  part,  this  reflects  an  initial  failure  to  appreciate  the 
importance  of  this  quantity,  which  will  now  be  discussed. 

To  consider  what  are  appropriate  values  of  Poisson’s  ratio  for  biological  tissues,  it  is 
helpful  to  consider  the  speed  of  longitudinal  and  shear  vibrations  in  tissue.  These  can  be  shown 
to  be: 


2  EQ-v) 

L  p(  l-2v)(l  +  v) 


(9) 


E 

2yO(l  +  v) 


(10) 


The  value  for  the  Young’s  modulus  E  of  normal  breast  tissue  is  on  the  order  of  10,000  Pa 
(van  Houten  et  al.,  [1999]  use  8000  Pa).  The  density  of  tissue  is  about  that  of  water,  so  p  is 
approximately  1000  kg/m3.  The  speed  of  longitudinal  vibrations  cL  is  just  the  sound  speed, 
which  is  about  1500  m/sec  in  tissue.  Using  Equation  (9),  we  can  readily  determine  that 
Poisson’s  ratio  must  be  very  close  to  0.5,  so  that  Equation  (9)  simplifies  to: 


E 

3yO(l-2v) 


(11) 


Equation  (11)  implies  that  v  -0.499999  if  the  longitudinal  sound  speed  is  1500  m/sec.  In  this 
case,  Equation  (10)  reduces  simply  to: 


(12) 


It  is  interesting  to  note  that  while  the  modulus  of  tissue  E  varies  over  a  reasonably  wide 
range,  the  longitudinal  sound  speed  is  found  to  be  quite  constant,  even  in  tumors  (the  constancy 
of  the  sound  speed  is  why  ultrasound  is  of  limited  use  for  detecting  cancer)  (Duck,  1990). 
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Equation  (11)  thus  implies  that  Poisson’s  ratio  must  vary  slightly  as  E  varies  to  keep  cL 
constant.  Arguably,  therefore,  a  more  appropriate  and  useful  parameterization  of  tissue  would 
use  Young’s  modulus  and  longitudinal  sound  speed  rather  than  Young’s  modulus  and  the 
Poisson’s  ratio. 

2.3.3  Numerical  Complications  Posed  by  Modeling  Poisson’s  Ratios  Near  0.5 

The  stiffness  matrix  KM  for  compressible  media  is  normally  developed  from  a  set  of 
constitutive  equations  that  relate  the  strains  e  to  stresses  a  (Smith  and  Griffiths,  1998).  In  the 
plane  strain  approximation,  the  associated  matrix  equation  is  written: 

a  =  De  (13) 


where  D  has  terms  of  the  form 


v 

T^v 


(14) 


As  v  approaches  0.5,  this  matrix,  and  thus  the  overall  stiffness  matrix  KM,  becomes  quite  ill- 
conditioned.  Special  techniques  must  be  used  to  obtain  a  converged  solution,  which  in  any  case 
can  be  difficult.  In  the  limit  where  the  Poisson’s  ratio  becomes  0.5,  the  material  becomes 
incompressible,  the  longitudinal  sound  speed  becomes  infinite,  and  a  different  formulation  of  the 
equations  of  motion  becomes  necessary.  In  this  formulation,  the  two  displacements  are  related 
by  the  incompressibility  condition  (conservation  of  volume),  and  for  this  reason  an  additional 
variable  must  be  introduced,  the  hydrostatic  pressure.  This  is  discussed  further  in  Section  2.4. 

Many  finite  element  codes  are  incapable  of  handling  the  purely  incompressible  case,  and 
often  have  convergence  difficulties  with  the  nearly-incompressible  case  as  well.  Numerical 
experiments  with  the  Creare-built  code  showed  poor  convergence  properties  as  Poisson’s  ratio 
was  increased  to  0.499  and  larger.  Presumably  for  the  same  reason,  much  of  the  MRE  work 
reported  in  the  literature  utilizes  values  of  0.45  to  0.49  (Van  Houten  et  al.,  2001;  Sinkus 
et  al.,  2000).  At  first  glance,  using  0.49  instead  of,  say,  0.49999,  appears  fairly  reasonable,  until 
one  calculates  the  sound  speeds.  At  a  value  of  0.49,  for  the  Young’s  modulus  and  density  cited 
earlier,  Equation  (9)  yields: 


cL~13m/scc  (15) 

This  is  far  less  than  the  actual  value  of  -1500  m/sec.  By  contrast,  the  shear  wave  speed,  given 
by  either  of  Equations  (10)  or  (12),  is  about  1.8  m/sec  and  is  insensitive  to  v . 

At  an  assumed  excitation  frequency  of  200  Hz,  in  the  range  typical  of  MRE  applications, 
the  corresponding  wavelengths  are  given  by  the  quotient  of  the  sound  speed  and  the  frequency: 

A*  ~  6.5 cm  (16) 

As  -0.9 cm  (17) 


*  Artificially  reduced  Poisson’s  ratio. 
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The  shear  wavelength  is  comparable  to  the  size  of  the  tumors.  As  mentioned  earlier,  achieving 
small  wavelengths  at  the  low  frequencies  resolvable  by  MRI  is  the  advantage  of  using  shear 
waves  in  the  first  place.  However,  the  calculated  longitudinal  wavelength  has  been  made  much 
shorter  than  the  actual  value  at  200  Hz  of  about  7  meters.  In  effect,  the  use  of  a  Poisson’s  ratio 
of  0.49  has  made  the  tissue  much  more  compressible  than  it  actually  is,  causing  the  wavelength 
of  longitudinal  waves  to  be  much  too  short. 

2.3.4  Errors  Caused  by  Artificial  Reductions  in  Poisson’s  Ratio 

As  mentioned  previously,  the  majority  of  the  published  work  on  MRE  has  involved 
artificial  reductions  in  Poisson’s  ratio,  just  as  in  the  initial  Creare  analyses,  presumably  to 
improve  the  numerical  performance  of  the  inverse  algorithm.  To  investigate  the  implications  of 
using  an  artificially  reduced  value  of  Poisson’s  ratio  on  the  calculated  tissue  behavior,  we 
performed  a  series  of  calculations  with  ANSYS.  This  code  was  used  rather  than  the 
corresponding  Creare  code  since  the  more  advanced  formulation  of  the  former  more  capably 
handles  nearly-incompressible  substances.  The  simulated  “breast”  in  these  calculations  is  a 
square  of  dimension  100  mm  x  100  mm,  excited  along  the  upper  edge  with  a  forced 
displacement  whose  amplitude  varies  sinusoidally  from  zero  at  the  left  and  right  comers  to  a 
maximum  value  of  1  mm  in  the  middle.  In  all  cases,  the  excitation  frequency  is  200  Hz,  the 
damping  constant/? is  .00016  (Gao,  1995),  and  the  density  p  is  1000  kg/m3.  The  nominal 

modulus  E  is  10000  Pa. 

Case  1:  Correct  Poisson’s  Ratio,  Longitudinal  Excitation 

The  first  case  involves  a  homogeneous  material,  and  a  stiff  inclusion  representing  a  tumor 
has  not  been  included.  The  Poisson’s  ratio  has  been  set  to  a  value  very  close  to  0.5  that  will 
yield  the  proper  longitudinal  sound  speed.  The  “breast”  is  excited  in  the  up-  and  down-  direction 
along  the  upper  edge  in  an  attempt  to  create  longitudinal  waves. 

A  “snapshot”  of  the  calculated  displacement  amplitudes  at  a  particular  time  is  shown  in 
Figure  3.  As  noted  above,  the  wavelength  of  a  200  Hz  longitudinal  wave  in  tissue  is  much  larger 
than  the  10  cm  maximum  dimension  of  the  breast.  Thus,  no  “waves”  are  seen,  just  positive 
amplitude  displacements  in  the  y-direction. 

Case  2:  Correct  Poisson’s  Ratio,  Shear  Excitation 

This  is  the  same  as  Case  1,  except  that  the  breast  is  excited  in  the  x-  (shear)  direction.  In 
this  case.  Figure  4  reveals  the  production  of  shear  waves  with  the  expected  wavelength  of  about 
0.9  cm  that  are  damped  fairly  quickly  with  depth.  Note  that  the  black  box  shown  in  the  figure 
has  no  significance  in  this  case  (the  box  denotes  the  location  of  the  inclusion  when  one  is  used  in 
Case  4). 

Case  3:  Artificially  Reduced  Poisson’s  Ratio,  Longitudinal  Excitation 

This  is  the  same  as  Case  1,  i.e.,  there  is  still  no  stiff  inclusion  and  the  breast  is  excited  in 
the  vertical  direction,  but  in  this  case  a  reduced  Poisson’s  ratio  of  0.49  is  used  as  in  some  of  the 
published  research.  Figure  5  exhibits  longitudinal  waves  with  a  wavelength  of  roughly  6.5  cm  as 
expected  from  Equation  (16).  The  comparison  of  Figure  5  to  the  correct  behavior  shown  in 
Figure  3  clearly  illustrates  the  adverse  effects  of  using  too  small  a  value  for  Poisson’s  ratio.  The 
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discrepancy  would  be  even  worse  if  still  smaller  values  of  v  were  used  (e.g.,  0.45),  as  is  seen  in 
the  literature. 

Case  4:  Correct  Poisson’s  Ratio,  Stiff  Inclusion,  Shear-Excitation 

This  case  is  similar  to  Case  2  in  that  it  employs  the  correct  Poisson’s  ratio.  Also,  in  this 
case  a  stiff  inclusion  is  provided  at  the  location  of  the  dark  box.  The  stiffness  (increased 
Young’s  modulus)  of  the  inclusion  leads  to  an  increase  in  sound  speed.  The  corresponding 
increase  in  shear  wavelength  across  the  inclusion  can  be  seen  in  Figure  6.  The  evident 
distortions  in  the  wave  patterns  provide  the  raw  data  necessary  for  an  inverse  method  to  estimate 
the  Young’s  modulus. 

As  demonstrated  by  the  results  shown  in  Figures  3  through  6,  the  assumed  Poisson’s  ratio 
exhibits  a  profound  effect  on  the  displacement  patterns.  While  the  use  of  a  reduced  value 
(e.g.,  0.49)  aids  numerical  convergence,  this  cannot  be  recommended  since  it  gives  rise  to 
relatively  low  longitudinal  wave  speeds,  correspondingly  short  wavelengths,  and  substantial 
distortion  of  the  calculated  displacements.  An  intriguing  question  is  whether  this  effect  may 
partly  explain  the  artifacts  reported  in  the  literature  when  in  vivo  data  and  that  from  phantoms 
are  used  to  reconstruct  the  stiffness.  To  date,  these  problems  have  been  attributed  to  the  use  of 
plane  strain  approximation  (rather  than  3-D)  or  to  the  assumption  of  linear  materials  models.  If 
more  detailed,  3-D  viscoelastic  models  are  needed,  this  will  severely  exacerbate  the  problems  of 
long-running  inverse  calculations  for  any  future  clinical  applications,  so  it  would  be  very  useful 
to  resolve  this  question. 

Since  the  longitudinal  wavelengths  at  the  frequencies  of  interest  are  much  larger  than  the 
dimensions  of  the  breast,  it  appears  justified  to  use  a  completely  incompressible  formulation  if 
this  will  improve  numerical  performance.  This  would  eliminate  the  numerical  problems 
associated  with  modeling  nearly-incompressible  behavior.  Such  an  approach  was  previously 
recommended  by  Skovoroda  et  al.  (1995).  An  approximate  incompressible  formulation  is 
developed  in  Section  2.4. 

2.3.5  Treatment  of  Damping 

Relatively  little  information  exists  in  the  literature  on  the  damping  of  low  frequency  shear 
waves  in  tissue.  What  information  does  exist  is  often  expressed  in  different  ways.  To  facilitate 
comparison  of  these  various  sources  of  information  and  to  support  the  use  of  different  codes  that 
require  the  damping  to  be  specified  in  various  ways,  the  following  expressions  were  developed, 
which  are  presented  here  without  proof: 

y 

a.  Damping  ratio  R:  - 

2  (op 

In  this  expression,  y  is  the  dimensional  damping  constant  (force  per  unit  volume  per  unit 
velocity),  (O  is  the  angular  frequency,  and  p  is  the  density.  According  to  Sinkus  et  al.  (2000), 
the  dimensional  damping  constant  of  tissue  is  typically  in  the  range  from  105  to  106  kg/m3-sec. 
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b.  Q  factor: 


Gao  (1995)  states  that  a  typical  value  of  Q  is  4,  which  is  consistent  with  the  R  values  provided  by 
Sinkus. 


c.  Attenuation  coefficient:  cc  = - 

2  csp 


In  this  expression,  cs  is  the  shear  wave  velocity, 
d.  Dimensionless  attenuation  coefficient:  c6 ls 

where  As  is  the  shear  wavelength.  Sinkus  et  al.  (2000)  state  that  typical  values  are  1-6  Np.  The 
relative  constancy  of  the  dimensionless  attenuation  coefficient  demonstrates  that  the  attenuation 
of  shear  waves  increases  strongly  with  frequency. 

2.4  Approximate  Approach  for  Modeling  An  Incompressible  Medium 

To  achieve  good  numerical  performance  without  introducing  modeling  artifacts  that  will 
result  from  reducing  Poisson’s  ratio,  an  appealing  solution  is  to  adopt  a  fully  incompressible 
formulation,  as  originally  recommended  by  Skovoroda  et  al.  (1995).  More  recently.  Bishop 
et  al.  (2000)  have  also  used  this  approximation  in  a  direct  inversion  scheme. 

One  can  motivate  this  approximation  in  several  ways.  Perhaps  the  simplest  way  is  to 
apply  Snell’s  law  to  an  incident  shear  wave  propagating  at  an  angle  6S  to  the  normal  of  an 

internal  tissue  interface: 


sin0L=—  sin^  (18) 

In  the  limit  where  the  ratio  of  the  longitudinal  and  shear  sound  speeds  is  very  large,  this 
expression  predicts  that  an  incident  shear  wave  will  not  produce  a  reflected  or  a  refracted 
longitudinal  wave  in  the  tissue  on  either  side  of  the  interface.  In  our  case,  the  speed  ratio  is 
roughly  1000,  so  neglecting  longitudinal  waves  seems  quite  justified. 

In  the  remainder  of  this  section,  we  will  develop  a  very  efficient,  but  approximate 
solution  for  the  incompressible  equations  of  motion. 

2.4.1  Derivation  of  Equations  of  Motion 

The  strain  and  stress  tensors,  respectively,  are  defined  by  Chenevert  et  al.  (1998): 
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and 

at=pSt+^Ee,  (20) 

where  the  i  and  j  indices  specify  coordinate  axes,  and  «,•  is  the  displacement  in  the  ith  direction. 
Specializing  to  the  2-D  plane  strain  approximation,  let  us  denote  by  u  the  displacements  in  the 
x-direction  and  use  v  to  refer  to  displacements  in  the  y-direction.  The  continuity  equation  in  2-D 
is  given  by: 


=  0  (21) 

dx  dy 


A  force  balance  on  a  small  volume  of  tissue  yields  the  equation  of  motion  in  the  x-direction: 


dii *  d  Mi 

,  a^a tmpir 


(22) 


where  y  is  the  (dimensional)  damping  coefficient  and  p  is  the  density. 
Equations  (19),  (20),  and  (21)  into  (22),  we  obtain  after  simplifying: 


£4* 


fa2M  aVi  du_  u 
dx2  dy2  j  Y  dt  P  dt2 


=  0 


Substituting 


(23) 


dp  +  1  ^ d2v  ^  d2v^  dv  d2' 
dy  3 


i J  v  v 

dx2  '  dy2frdt~P~di2  = 


(24) 


Note  that  in  developing  Equations  23  and  24,  we  have  neglected  terms  involving  spatial 
derivatives  of  the  Young’s  modulus.  The  validity  of  this  key  simplification  will  be  discussed 
below. 

We  can  eliminate  the  pressure  by  subtracting  from  the  derivative  of  Equation  23  with 
respect  to  y  the  derivative  of  Equation  24  with  respect  to  x.  Next,  in  analogy  with  a  very  similar 
concept  introduced  in  incompressible  fluid  mechanics,  we  define  the  “vorticity”  by: 


£ 


du  dv 
dy  dx 


(25) 


After  manipulation,  we  obtain: 


3  [dx2  dy2)  Y  dt  Pdt2 


=  0 


(26) 


If  necessary,  the  displacements  can  be  recovered  by  introducing  a  function  y/  defined  as 
follows: 
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d\f/ 

U~~dy 

dy/ 

V"  17 

The  quantity  yf  is  called  the  stream  function  in  hydrodynamics,  and  its  use  automatically  ensures 
that  the  continuity  Equation  (21)  is  always  satisfied.  Using  the  definition  Equation  (25),  the 
stream  function  is  given  by: 
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(28) 


avaV 

*  dx2  dy2 


(29) 


Thus,  one  can  solve  the  incompressible  equations  of  motion  using  the  following 
sequential  procedure: 

a.  Calculate  the  vorticity  using  Equation  (26) 

b.  Calculate  the  stream  function  using  Equation  (29) 

c.  Differentiate  the  stream  function  to  yield  the  displacements  via  Equations  (27)  and  (28). 

We  are  not  especially  interested  in  the  displacements  per  se,  but  in  what  they  can  tell  us 
about  the  Young’s  modulus.  This  being  the  case,  a  potentially  more  efficient  procedure  would 
be  to  base  the  inverse  problem  not  on  the  displacements  themselves,  but  on  their  vorticity.  This 
eliminates  the  need  to  solve  Equations  (29),  (27)  and  (28).  To  do  this  one  would  first  obtain 
estimates  for  the  vorticity  distribution  by  differentiating  the  MRI-supplied  displacement  field, 
possibly  smoothing  the  results  afterward  to  eliminate  high  frequency  measurement  noise.  Then 
one  would  determine  a  set  of  Young’s  moduli  that  provide  the  best  fit  between  the 
experimentally  derived  vorticities  and  those  calculated  by  Equation  (26).  This  would  be 
accomplished  by  inputting  to  Equation  (26)  the  experimentally  measured  vorticity  values  at  the 
boundaries  of  the  breast,  and  iterating  the  assumed  values  of  the  Young’s  moduli  until  a  good 
match  is  obtained  to  the  interior  vorticities.  The  details  of  how  such  an  inverse  problem  can  be 
solved  are  discussed  in  Section  2.5. 

We  note  in  passing  that  a  somewhat  more  general  derivation  of  Equation  (26)  in  the 
absence  of  damping  is  provided  by  Landau  and  Lifshitz  (1986). 

2.4.2  Preliminary  Assessment  of  the  Vorticity  Approach 

Let  n  be  the  number  of  grid  points  used  for  discretizing  the  motion  of  the  breast  tissue  in 
2-D.  The  fundamental  equations  of  motion,  Equations  (21)  involve  3n  variables,  the  pressure 
and  the  displacements  along  both  axes.  These  can  be  obtained  by  solving  the  3n  Equations  (21) 
and  (22)  (it  is  also  possible  to  reduce  the  number  of  equations  to  2n  by  formally  eliminating 
pressure,  as  described  by  Bishop  et  al.,  2000).  The  prime  advantage  of  the  vorticity  procedure  is 
that  the  displacement  field  can  be  characterized  by  solving  only  n  Equations  (26).  The  time 
required  to  obtain  a  solution  goes  roughly  as  the  cube  of  the  number  of  equations,  so  this 
represents  a  potential  saving  in  computer  time  of  a  factor  between  eight  and  twenty-seven,  a  very 
significant  figure.  A  somewhat  more  subtle  benefit  of  the  vorticity  formulation  is  that  the 
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underlying  equations  are  symmetric,  i.e.,  are  unchanged  if  the  indices  i  and  j  are  exchanged.  As 
mentioned  earlier,  efficient  solution  techniques  are  available  in  this  case  for  solving  the 
associated  matrix  equation.  By  contrast  the  3n  set  of  Equations  (21)  and  (22)  are  not  symmetric 
and  a  less  efficient  solution  technique  must  be  used. 

However,  two  potential  disadvantages  are  associated  with  the  vorticity  approach.  First, 
the  derivatives  of  the  MRI-measured  displacement  data  are  used,  both  as  boundary  conditions 
and  to  serve  as  the  basis  of  comparison  (or  “figure-of-merit”)  used  to  solve  the  inverse  problem. 
Differentiation  of  the  data  will  likely  increase  the  sensitivity  to  noise  in  this  data.  Second,  as 
mentioned  above,  terms  involving  the  spatial  derivatives  of  the  Young’s  modulus  were  neglected 
in  the  derivation  of  the  Equations  (26).  As  discussed  below,  these  terms  may  not  always  be 
negligible  if  very  steep  gradients  in  tissue  stiffness  are  encountered. 

2.4.3  Implementation  of  the  Vorticity  Approach 

To  facilitate  experimentation  with  the  vorticity  approach,  a  finite  element  formulation  of 
Equations  (26)  was  developed  in  Fortran  90  using  the  same  techniques  employed  for  the 
compressible  displacement-based  approach  discussed  in  Section  4.  As  before,  the  two- 
dimensional  problem  domain  was  taken  to  be  square  (rather  than  a  circle)  for  simplicity,  and 
9-node  quadrilateral  elements  were  used. 

* 

To  obtain  a  solution,  a  calculation  was  initially  made  with  ANSYS,  which  in  effect  took 
the  place  of  an  MRI  machine.  These  calculations  assumed  that  the  upper  boundary  of  the  square 
was  excited  in  the  transverse  direction  by  a  driving  function  whose  amplitude  varied  sinusoidally 
in  the  x-direction.  The  other  boundaries  were  assumed  to  be  free.  The  resulting  displacements 
were  differentiated  numerically  to  calculate  vorticities.  The  values  of  these  ANSYS-derived 
vorticities  around  the  four  edges  of  the  square  were  input  as  boundary  conditions  to  the  Creare 
code.  If  MRI  data  was  used,  the  same  procedure  would  be  followed,  except  that  filtering  might 
be  necessary  to  reduce  the  sensitivity  to  noise.  Utilizing  displacement  data  around  the  edges 
allows  the  simulation  to  at  least  partially  reflect  the  edge  restraint  provided  by  adjacent  tissue 
layers. 


The  Creare  vorticity  code  executes  very  rapidly  for  a  given  assumed  distribution  in 
Young’s  modulus.  For  example,  a  solution  can  be  obtained  on  a  100  x  100  grid  in  only  about 
three  seconds  using  a  1  GHz  personal  computer.  Such  a  grid  allows  a  planar  slice  of  the 
simulated  breast  to  be  discretized  relatively  finely,  as  squares  of  about  0.1  x  0.1  mm.  This  is 
considered  adequate  for  even  relatively  high  frequencies,  and  is,  in  fact,  considerably  finer  than 
many  of  the  simulations  reported  in  the  literature. 

2.4.4  Comparison  of  the  Vorticity  Approach  to  a  Conventional  Displacement-Based 
Scheme 

As  shown  in  Figures  7  and  8,  very  good  agreement  was  obtained  between  the  ANSYS 
results  and  those  obtained  with  the  Creare  code  for  cases  with  a  spatially  constant  Young’s 
modulus.  It  was  found  that  the  agreement  was  improved  if  the  modulus  used  in  the  latter  was 
increased  by  about  seven  percent  over  the  value  assumed  in  ANSYS;  this  is  attributed  to 
relatively  minor  differences  in  numerical  technique. 
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However,  larger  differences  were  observed  when  very  strong  discontinuities  in  Young’s 
modulus  were  introduced  to  simulate  the  presence  of  a  tumor.  In  the  calculations  discussed 
below,  the  tumor  was  represented  by  an  area  that  has  a  Young’s  modulus  that  is  ten  times  that  of 
healthy  tissue  (100  kPa  versus  the  nominal  stiffness  of  10  kPa).  The  location  of  the  tumor  was  in 
the  upper  right-hand  comer  of  the  problem  domain,  i.e.  at: 

0.02 m  <x<  .03m 

0.02m  <  y  <  .03 m 

Figures  9  and  10  show  the  vorticities  calculated  by  ANSYS  and  the  Creare  code  for  this 
case.  While  these  do  not  appear  too  different,  differences  can  be  seen  in  the  detailed  behavior 
around  the  inclusion.  To  see  this  even  more  clearly,  consider  a  similar  case  where  the  excitation 
frequency  is  reduced  from  200  Hz  to  50  Hz.  The  overall  vorticities,  Figures  11  and  12,  again 
appear  qualitatively  similar.  However,  Figure  13  shows  the  vorticity  values  along  a  vertical  path 
through  the  simulated  tumor.  Note  that  a  sharp  discontinuity  occurs  at  the  bottom  edge  of  the 
tumor  where  a  sudden,  large  change  in  Young’s  modulus  is  experienced.  A  smaller  discontinuity 
occurs  at  the  top.  These  jumps  in  vorticity  can  be  shown  to  result  from  terms  involving  the 
second  spatial  derivatives  of  the  Young’s  modulus  that  were  neglected  in  the  derivation  of 
Equations  (26).  While  qualitatively  similar,  the  vorticities  calculated  by  the  Creare  code  along 
the  same  path  have  no  such  discontinuities.  The  differences  in  vorticity  near  the  excited  edge 
(y=.05  cm)  are  believed  to  be  caused  by  differences  in  reflections  at  the  tumor. 

For  a  similar  problem  in  which  the  (again  discontinuous)  change  in  Young’s  modulus  is 
assumed  to  be  “only”  a  factor  of  2  and  the  frequency  is  restored  to  the  nominal  value  of  200  Hz, 
the  “vorticity  jump”  calculated  by  ANSYS  is  much  less  pronounced,  as  shown  on  an  expanded 
scale  in  Figure  14.  Also,  the  results  of  the  two  codes  in  this  case  agree  fairly  well.  The  vorticity 
jump  would  be  smaller  still  if  the  tissue  properties  were  allowed  to  vary  more  smoothly  than  the 
artificial  step  changes  we  have  assumed.  It  remains  to  be  seen  whether  vorticity  jumps  in  actual 
tissue  are  so  large  as  to  compromise  our  ability  to  estimate  Young’s  modulus  using  this  method. 
If  the  jumps  are  indeed  that  dramatic,  we  can  speculate  that  the  vorticity  concept  could  provide  a 
means  to  process  raw  displacement  data  and  detect  sudden  changes  in  Young’s  modulus  without 
solving  an  inverse  problem  at  all. 

2.4.5  Discussion 

An  extremely  efficient  approximate  solution  technique  for  the  tissue  equations  of  motion 
was  developed.  Using  a  relatively  modest  1  GHz  personal  computer,  a  solution  can  be  obtained 
on  a  100  x  100  planar  grid  in  only  about  three  seconds.  Further  increases  in  execution  speed  are 
possible. 

Two  key  assumptions  contribute  to  the  efficiency  of  this  solution.  First,  the  tissue  is 
assumed  to  be  incompressible.  This  implicitly  neglects  the  existence  of  longitudinal  waves, 
which  is  justified  on  the  basis  of  their  relatively  high  propagation  velocities  and  long 
wavelengths.  Second,  the  Young’s  modulus  is  assumed  to  vary  sufficiently  so  slowly  that  its 
higher  order  spatial  derivatives  can  be  neglected. 
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Further  work  is  necessary  to  validate  whether  the  vorticity  approximation  can  be  used 
when  analyzing  actual  tissue  displacements.  If  the  vorticities  calculated  from  MRI  displacement 
data  exhibit  large,  discontinuous  jumps,  this  may  invalidate  the  use  of  vorticity  for  quantitatively 
modeling  tissue  response.  However,  in  this  case,  the  vorticity  concept  could  still  be  useful  since 
the  vorticities  calculated  from  displacement  data  may  provide  a  direct  and  even  more  efficient 
means  to  locate  regions  where  tissue  stiffness  is  rapidly  varying. 

2.5  Solution  of  the  Inverse  Problem 
2.5.1  Overview 

A  conventional  scheme  for  solving  the  inverse  problem  is  to  minimize  a  figure  of  merit 
of  the  form: 


^  =  Z  ~  ) 2  +  o.  -  VD2  (3°) 

In  Equation  (30),  u  and  v  denote  the  x-  and  y-direction  displacements  at  each  of  the  n 
nodes  of  a  plane  strain  displacement  model  at  some  specified  time  in  their  oscillatory  cycle.  The 
superscript  “°”  denotes  displacements  obtained  from  measured  MRI  data  (which  can  come  from 
an  actual  MRI  or  from  a  finite  element  simulation).  The  lack  of  a  superscript  denotes  code- 
calculated  displacements  based  on  the  current  estimate  for  Young’s  modulus  at  each  node.  The 
same  expression  would  be  used  for  a  3-D  analysis  or  one  based  on  vorticity,  the  only  difference 
being  the  number  of  terms  included  in  the  sum. 

The  Dartmouth  group  determines  Young’s  modulus  by  minimizing  J  using  a  well-known 
technique  known  as  the  Levenberg-Marquardt  algorithm  (Press  et  al.,  1986).  This  requires  the 
estimation  of  the  Hessian  matrix  associated  with  J,  which  in  turn  requires  the  calculation  of 
terms  such  as: 


dv 

dE 


i=l,...n;  k=l,...p 


(31) 


where  E *  denotes  the  Young’s  modulus  of  the  kth  element  in  the  finite  element  model.  If  there 
are  n  nodes  and  p  elements  in  the  model,  the  calculation  of  the  nodal  displacements  requires  a 
time  proportional  to  roughly  n3.  However,  the  calculation  of  the  terms  (31)  requires  a  time  that 
is  p  times  longer  than  this,  which  for  low  order  elements  (for  which  p  is  roughly  equal  to  n)  is  on 
the  order  of  n4.  This  is  a  significant  disadvantage  as  the  number  of  elements  and  nodes  grows 
unless  the  calculation  of  these  quantities  makes  convergence  to  the  proper  Young’s  modulus 
distribution  much  more  rapid. 

Since  the  time  required  to  calculate  the  terms  in  the  approximate  Hessian  matrix  is  very 
long,  we  instead  investigated  the  use  of  a  quasi-Newton  method,  specifically  the  Boyden- 
Fletcher-Goldfarb-Shanno  (BFGS)  algorithm  (Press  et  al.,  1986).  While  this  method  is  likely  to 
require  more  iterations  to  converge  than  Levenberg-Marquardt,  BFGS  requires  the  calculation  of 
only  the  gradient  of  the  figure-of-merit,  i.e.,  the  terms: 
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To  facilitate  use  of  this  technique,  we  developed  an  algorithm  that  very  efficiently  calculates  the 
terms  of  the  gradient. 

2.5.2  Calculation  of  the  Gradient  of  the  Figure-of-Merit 

As  discussed  previously,  given  a  set  of  measured  (or  simulated)  displacement  or  vorticity 
data,  one  adjusts  the  assumed  Young’s  moduli  to  make  the  calculated  values  agree  with  the 
measured  quantities.  This  can  be  formulated  as  an  optimization  problem,  in  which  a  measure  of 
the  disagreement,  Equation  (30),  is  minimized. 

Most  efficient  and  robust  techniques  for  minimizing  figure-of-merit  functions  of  the  form 
of  Equation  (30)  require  information  on  the  change  in  the  predicted  displacements  caused  by  a 
small  change  in  input  parameters.  For  example,  the  Levenberg-Marquardt  technique  requires 
extensive  information  that  is  time-consuming  to  calculate,  i.e.,  the  terms  shown  in  Equation  (31). 
The  BFGS  technique  requires  less  gradient  information,  but  it  is  still  critical  to  develop  this 
information  in  an  efficient  manner.  For  example,  the  crudest  approach  to  calculate  the  gradient 
terms  shown  in  Equation  (32)  would  involve  making  a  small  change  in  a  single  modulus  E*  and 
then  recalculating  the  entire  displacement  solution.  Repeating  this  process  p  times  to  calculate 
all  such  terms  would  again  require  a  time  on  the  order  of  p  times  n3,  or  roughly  n  .  In  other 
words,  this  brute-force  process  would  require  the  same  length  of  time  as  the  Levenburg- 
Marquardt  method. 

Note  that  the  brute-force  technique  calculates  the  change  in  each  displacement  or 
vorticity  caused  by  a  change  in  the  input  Young’s  moduli.  Since  we  then  sum  over  these 
displacements  to  calculate  the  figure-of-merit  J,  we  see  that  the  brute-force  technique  calculates 
much  more  information  than  we  really  need  for  BFGS.  A  better  scheme  is  to  use  an  adjoint 
technique  that  efficiently  calculates  only  the  needed  information  (Kenton,  1981).  If  we  have  a 
set  of  algebraic  equations  of  the  form: 


0  =  g,(x) 


(33) 


together  with  a  single  figure-of-merit  J(x),  one  defines  an  adjoint  x+  associated  with  each 
variable  x  by  constructing  the  following  equations: 


o=Eir*;  +f- 

k  OXj 


(34) 


After  the  terms  x+  have  been  calculated  using  Equation  (34),  the  terms  of  Equation  (32)  can  be 
obtained  very  simply  by  calculating  the  following  sums: 
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In  our  case,  Equations  (33)  are  linear,  and  the  matrix  terms  in  Equation  (34)  put  into  the  standard 
form  of  Equation  (6)  are  of  the  form: 

Aik  =  =  KMki  +  ifiaiKM ki  -  tu2MMki  (36) 

dxi 

When  the  stiffness  and  mass  matrices  KM  and  MM  are  symmetric,  and  when  this  symmetry  is 
not  destroyed  by  using  an  injudicious  formulation  of  the  boundary  conditions,  the  matrix  A'  used 
to  solve  for  the  adjoints  is  the  same  as  the  matrix  A  used  to  solve  for  the  displacements  or  the 
vorticities. 

This  symmetry  offers  the  potential  for  a  still  larger  savings  in  computation  time,  if 
sufficient  memory  is  available  to  use  direct  sparse  matrix  techniques  rather  than  an  iterative 
scheme  such  as  PCG.  If  we  use  such  a  direct  matrix  solution  technique,  we  will  manipulate  A  in 
some  fashion  (e.g.,  decompose  it  into  upper-  and  lower-diagonal  matrices)  that  will  then  allow 
the  matrix  equation  to  be  readily  solved.  The  time-consuming  part  of  the  solution  process  is  the 
matrix  manipulation  step,  which  takes  a  time  on  the  order  of  the  cube  of  the  rank,  i.e.,  in  this  case 
n3.  The  subsequent  “back-substitution”  process  needed  to  obtain  the  solution  *  for  a  given 
b  takes,  by  contrast,  only  a  time  on  the  order  of  n2.  Thus,  direct  methods  are  very  appealing 
when  used  with  the  BFGS  algorithm:  after  the  matrix  A  has  been  manipulated  to  obtain  a 
solution  of  the  forward  problem,  the  adjoints  can  be  obtained  with  the  same  matrix 
decomposition  in  essentially  negligible  additional  time.  The  gradient  terms  are  then  obtained 
trivially  using  Equation  (35). 

2.5.3  Demonstration  Calculation  of  Determining  Young’s  Modulus  From  Data 

To  experiment  with  the  inverse  algorithm,  both  the  compressible,  displacement-based 
finite  element  code  and  the  incompressible  vorticity  code  were  augmented  to  calculate  the 
gradient  terms  using  the  adjoint  method.  After  doing  this,  it  was  first  confirmed  that  the  change 
in  figure-of-merit  caused  by  a  change  in  one  or  more  of  the  Young’s  moduli  calculated  using 
brute-force  techniques  was  the  same  as  predicted  by  the  adjoint  method.  The  gradient 
information  was  then  used  to  minimize  the  figure-of-merit  Equation  (30).  We  tried  several 
algorithms  to  achieve  this,  and  the  best  results  were  obtained  using  a  commercially-available 
implementation  of  the  BFGS  algorithm  (IMSL,  1994).  Figure  15  shows  the  relative  rate  at 
which  the  figure-of-merit  converged  in  a  typical  problem  (a  value  of  zero  for  the  figure-of-merit 
implies  perfect  convergence). 

Figure  16  shows  the  Young’s  modulus  obtained  from  computer-generated  MRI 
displacement  vorticity  “data.”  Computer-generated  data  were  obtained  by  merely  running  the 
finite  element  code  on  a  100  x  100  grid  for  a  specified  Young’s  modulus  distribution,  in  this  case 
a  uniform  distribution  with  a  value  of  10  kPa  except  for  a  square  “tumor”  of  value  100  kPa  in  the 
upper-right-hand  comer.  The  minimization  procedure  then  begins  with  a  uniform  (and  incorrect) 
Young’s  modulus  assumption,  and  iterates  until  the  desired  accuracy  is  obtained.  As  can  be  seen 
in  Figure  16,  the  algorithm  does  an  excellent  job  inferring  Young’s  modulus  from  the  computer¬ 
generated  data.  This  simulation  required  approximately  10  minutes  of  computer  time  on  a 
1000  MHz  personal  computer.  While  several  possibilities  exist  for  further  increasing  the  speed 
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of  this  procedure,  the  relatively  short  time  required  to  obtain  good  convergence  is  a  direct 
consequence  of: 

a.  Use  of  the  vorticity  formulation,  which  results  in  a  greatly  reduced  number  of 
equations  for  a  given  discretization. 

b.  The  use  of  a  boundary  condition  formulation  that  results  in  a  symmetric  A  matrix. 

c.  Exploiting  the  adjoint  formulation  to  very  efficiently  calculate  the  gradient  of  the 
figure  of  merit. 

As  shown  in  Figure  17,  similar  success  was  also  obtained  using  a  displacement-based 
formulation  (with  an  artificially  reduced  Poisson’s  ratio).  Should  the  vorticity  scheme  prove  too 
inaccurate  when  applied  to  actual  MRI-measured  tissue  displacements,  this  success  indicates  that 
an  incompressible  version  of  the  displacement  model  such  as  that  used  by  Bishop  et  al.  (2000) 
should  converge  just  as  well  as  the  vorticity  scheme,  but  in  a  longer  time.  Nevertheless,  even  in 
this  case  the  efficiency  advantages  of  using  adjoints  should  permit  convergence  in  a  much 
shorter  period  than  is  required  by  other  techniques. 

Acknowledging  the  difficulties  implicit  in  matching  the  vorticity  model  to  ANSYS 
results  when  sharp  discontinuities  in  assumed  stiffness  are  present,  a  final  inverse  problem  was 
run  in  which  the  ANSYS  results  were  input  to  the  vorticity  approach.  The  reconstructed 
Young’s  modulus  distribution  is  shown  in  Figure  18  (note  that  the  figure  actually  shows 
“smoothed”  data,  in  which  each  modulus  value  was  averaged  with  its  nearest  neighbors  to 
suppress  high  frequency  oscillations  that  occur  near  the  vorticity  jumps).  The  results  are 
considered  surprisingly  good,  and  the  tumor  can  clearly  be  discerned.  This  suggests  that  even 
with  its  simplifying  assumptions,  the  vorticity  approach  may  still  yield  acceptable  results  when 
applied  to  real  data. 

One  important  detail  of  the  inverse  algorithm  should  be  explained.  Initial  numerical 
experiments  resulted  in  a  much  poorer  reconstruction  of  the  Young’s  modulus  than  is  shown  in 
Figures  16  and  17.  Inevitably,  we  found  that  the  method  converged  on  a  solution  that  involved 
non-physical  (negative)  Young’s  moduli  adjacent  to  the  “tumor.”  Further,  the  calculated 
Young’s  moduli  inside  the  tumor  boundary  compensated  for  this  and  were  too  large.  Other 
numerical  artifacts  involving  negative  Young’s  moduli  also  sometimes  appeared  near  the 
boundary  of  the  simulated  breast.  To  solve  these  problems  without  instituting  the  full 
complication  of  a  constrained  minimization  approach,  the  Young’s  moduli  were  mapped  to  a 
new  variable  a ,  and  convergence  was  obtained  on  a.  The  mapping  used: 

(37) 

which  guarantees  that  the  Young’s  modulus  will  be  larger  than  some  physically-based  minimum 
value  Emin.  This  modification  to  the  algorithm  proved  very  beneficial. 

In  these  demonstration  calculations,  the  damping  was  assumed  known.  When  applied  to 
real  data,  the  damping  may  well  vary  across  a  cross  section  of  the  breast,  although  there  is  little 
in  the  literature  to  imply  how  much  variation  should  be  expected.  In  any  event,  simultaneously 
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solving  for  unknown  damping  coefficients  along  with  the  Young’s  moduli  should  not  lead  to  any 
additional  complications.  In  particular,  the  gradient  of  the  figure-of-merit  with  respect  to  the 
damping  is  obtained  from  an  expression  of  the  form  of  Equation  (35),  using  the  same  values  for 
the  adjoint  functions  as  are  used  to  calculate  the  gradient  with  respect  to  the  Young’s  moduli. 


24 


TM-2201 


<§reare 

3  KEY  ACCOMPLISHMENTS 

•  An  adjoint  method  developed  on  this  project  provides  an  extremely  efficient 
technique  for  calculating  the  gradient  of  the  goodness-of-fit  measure  that  is  the  basis 
of  the  inverse  algorithm. 

•  Combining  the  adjoint  method  with  the  BFGS  algorithm  provides  an  efficient  means 
for  calculating  tissue  stiffness  from  MRI  displacement  data. 

•  Much  of  the  experience  with  Magnetic  Resonance  Elastography  previously  reported 
in  the  literature  has  utilized  an  incorrect  value  of  Poisson’s  ratio. 

•  Use  of  an  incorrect  Poisson’s  ratio  leads  to  a  drastic  underprediction  of  the 
longitudinal  acoustic  wavelength,  which  in  turn  can  be  expected  to  cause  artifacts  in 
the  reconstructed  hardness  values. 

•  An  incompressible  tissue  model  developed  under  this  project  provides  an  extremely 
efficient  method  for  calculating  tissue  displacements. 

•  Demonstration  of  the  speed  and  accuracy  of  our  incompressible  tissue  model  linked 
to  our  inverse  problem  solution  method  using  2-D  computer-generated  data. 
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4  REPORTABLE  OUTCOMES 

•  Era  of  Hope  Poster  Presentation  for  Department  of  Defense  Breast  Cancer  Research 
Program  Meeting,  September  25,  2002,  Orlando,  FL.  (See  Abstract  at  end  of 
Appendices  section.) 

•  Grant  entitled  “MR  Elastography  for  Breast  Cancer  Detection,”  submitted  to 
NIH/SBIR  July  25,  2002. 
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5  CONCLUSIONS 

In  this  report,  we  have  developed  an  efficient  method  for  determining  Young’s  modulus 
from  tissue  displacement  data.  The  method  is  based  on  minimizing  a  figure-of-merit  that 
indicates  how  closely  the  experimental  measurements  match  a  model  of  the  tissue  response.  One 
hallmark  of  this  method  is  the  use  of  an  adjoint  technique  for  calculating  the  gradient  of  the 
figure-of-merit  with  respect  to  the  unknown  quantities  (i.e.,  Young’s  modulus  and,  perhaps, 
damping)  in  a  time  that  is  negligible  compared  to  the  solution  of  the  tissue  equations  of  motion. 

Use  of  artificially  reduced  values  for  the  Poisson’s  ratio,  as  is  often  done,  was  found  to 
introduce  severe  artifacts  into  the  calculated  displacements.  A  better  approach  is  to  assume 
completely  incompressible  behavior. 

Very  fast  convergence  of  the  calculated  Young’s  modulus  distribution  was  obtained  with 
simulated  displacement  data  by  using  a  “vorticity”  scheme  that  models  2-D  incompressible 
behavior  using  a  greatly  reduced  number  of  equations.  Further  work  is  needed  to  determine 
whether  the  approximations  inherent  in  this  approach  will  be  sufficiently  accurate  when  applied 
to  actual  in  vivo  data. 
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Figure  1 .  Comparison  of  Finite  Element  Code  Predictions  of  the  Logarithm  of  the  Displacement 
Amplitude  to  the  Logarithm  of  the  Displacement  Amplitude  for  the  Exact  Solution.  These  data  show  that 
our  finite  element  code  closely  matches  the  exact  solution  for  the  analyzed  configuration. 


Figure  2.  Instantaneous  Displacements  Versus  Depth  for  Various  Meshing  Schemes.  These  data  show 
the  effect  of  model  order  on  the  match  to  the  exact  solution.  For  this  model,  3600  quadrilateral  elements 
or  20,000  triangular  elements  yield  good  model  accuracy. 
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Figure  3.  Calculated  Longitudinal  Displacements  for  Longitudinal  Excitation  With  Correct  Material 
Properties  Using  ANSYS. 
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Figure  4.  Calculated  Shear  Displacements  for  Shear  Excitation  With  Correct  Material  Properties  Using 
ANSYS. 
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Figures.  Calculated  Longitudinal  Displacement  for  Longitudinal  Excitation  With  Artificially  Lowered 
Poisson’s  Ratio  Using  ANSYS. 
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Figure  6.  Calculated  Shear  Displacements  for  Shear  Excitation  With  Correct  Material  Properties  Using 
ANSYS  with  Stiff  Inclusion. 
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Figure  7.  Vorticities  Obtained  by  Differentiating  ANSYS-Calculated  Displacements  for  a  Homogenous 
Material  (constant  stiffness)  Excited  at  200  Hz. 
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Figure  9.  Vorticities  Obtained  by  Differentiating  ANSYS-Calculated  Displacements  for  the  Case  of  a 
Simulated  Tumor  Imaged  at  200  Hz. 


Figure  10.  Vorticities  Calculated  Using  Creare  Finite  Element  Code  for  the  Case  of  a  Simulated  Tumor 
Imaged  at  200  Hz. 
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Figure  11.  Vorticities  Obtained  by  Differentiating  ANSYS-Calculated  Displacements  for  Case  With  a 
Simulated  Tumor  Imaged  at  50  Hz. 


Figure  12.  Vorticities  Calculated  Using  Creare  Finite  Element  Code  for  Case  of  a  Simulated  Tumor 
Imaged  at  50  Hz. 
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Figure  13.  Comparison  of  Creare  and  ANSYS  Vorticity  Values  for  a  Case  With  a  Sharp  Discontinuity  in 
Young’s  Modulus  at  Tumor  Boundaries  of  y=.02  and  .03  cm. 


0.010  0.015  0.020  0.025  0.030  0.035  0.040 


-  Y,  cm 

•  ANSYS 

▲  Creare  Vorticity  Code 

Figure  14.  Vorticities  for  a  Case  Involving  a  Tumor  With  Only  Twice  the  Nominal  Stiffness,  Imaged  at 
200  Hz;  Little  Discontinuity  in  ANSYS-Calculated  Vorticity  Is  Evident  at  the  Tumor  Boundaries  of  y=.02 
and  y=.03  cm,  and  the  Two  Codes’  Results  Agree  Fairly  Well. 
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Figure  17.  Young’s  Modulus  Calculated  from  Computer-Generated  MRI  Displacement  Data  Using  the 
Creare  Displacement-Based  Code  and  Inverse  Algorithm.  The  stiffness  (relative  modulus  10)  of  the 
simulated  tumor  in  the  upper-right-hand  corner  can  be  easily  discerned. 
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Figure  18.  Reconstructed  (averaged)  Young’s  Modulus  Using  Displacement  Data  From  ANSYS 
Calculation  as  Input  to  Vorticity  Model.  Despite  the  differences  between  the  two  approaches,  the 
presence  and  location  of  the  tumor  are  readily  apparent 
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One  physical  characteristic  that  clearly  distinguishes  healthy  from  cancerous  tissue  is  mechanical 
stiffness  (hardness).  The  Young’s  modulus  of  breast  tumors  can  be  one  to  two  orders  of 
magnitude  larger  than  that  of  normal  tissue.  For  this  reason,  palpation  has  long  been  used  by 
physicians  and  patients  for  early  screening.  However,  manual  palpation  is  not  very  effective  for 
tumors  lying  deep  within  the  breast  and  is  not  quantitative.  For  these  reasons,  over  the  past 
decade,  numerous  efforts  have  been  underway  to  combine  external  mechanical  stimulation  and 
various  imaging  techniques  to  quantitatively  measure  the  Young’s  modulus  of  tissue  throughout 
both  the  breast  and  the  prostate.  This  has  been  termed  “palpation  at  a  distance.”  Because  of  its 
inherent  sensitivity,  much  current  research  is  focused  on  using  magnetic  resonance  imaging 
(MRI),  and  the  overall  process  is  termed  magnetic  resonance  elastography  (MRE). 

One  of  the  most  challenging  technical  aspects  of  MRE  is  the  solution  of  the  “inverse  problem,” 
i.e.,  quantitatively  determining  Young’s  modulus  from  MRI-measured  tissue  displacement  data. 
Creare  is  developing  analytical  solution  techniques  to  improve  the  efficiency  and  robustness  of 
the  inverse  problem  solution.  We  are  also  investigating  improvements  in  the  physical  modeling 
on  which  the  inverse  method  is  based,  primarily  in  the  representation  of  tissue  compressibility. 

During  this  project,  we  have  shown:  (1)  an  adjoint  method  developed  on  this  project  provides  an 
extremely  efficient  technique  for  calculating  the  gradient  the  goodness-of-fit  measure  that  is  the 
basis  of  the  inverse  algorithm;  (2)  combining  the  adjoint  method  with  a  quasi-Newton  method, 
specifically  the  Boyden-Fletcher-Goldfarb-Shanno  (BFGS)  algorithm,  provides  an  efficient 
means  for  calculating  tissue  stiffness  from  MRI  displacement  data;  (3)  much  of  the  experience 
with  MRE  previously  reported  in  the  literature  has  utilized  an  incorrect  value  of  Poisson’s  ratio; 
and  (4)  use  of  an  incorrect  Poisson’s  ratio  leads  to  a  drastic  underprediction  of  the  longitudinal 
acoustic  wavelength,  which  in  turn  can  be  expected  to  cause  artifacts  in  the  reconstructed 
hardness  values. 
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